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Abstract 

Background: The hypothalamic-pituitary-adrenal axis (HPA axis) is a major part of the 
neuroendocrine system responsible for the regulation of the response to physical or 
mental stress and for the control of the synthesis of the stress hormone Cortisol. 
Dysfunctions of the HPA axis characterized by either low (hypocortisolism) or increased 
(hypercortisolism) Cortisol levels are implicated in various pathological conditions. Their 
understanding and therapeutic correction may be supported by mathematical 
modeling and simulation of the HPA axis. 

Methods: Mass action and Michaelis Menten enzyme kinetics were used to provide a 
mechanistic description of the feedback mechanisms within the pituitary gland cells by 
which Cortisol inhibits its own production. A separation of the nucleus from the 
cytoplasm by compartments enabled a differentiation between slow genomic and fast 
non-genomic processes. The model in parts was trained against time resolved ACTH 
stress response data from an in vitro cell culture of murine AtT-20 pituitary tumor cells 
and analyzed by bifurcation discovery tools. 

Results: A recently found pituitary gland cell membrane receptor that mediates rapid 
non-genomic actions of glucocorticoids has been incorporated into our model of the 
HPA axis. As a consequence of the distinction between genomic and non-genomic 
feedback processes our model possesses an extended dynamic repertoire in 
comparison to existing HPA models. In particular, our model exhibits limit cycle 
oscillations and bistable behavior associated to hypocortisolism but also features a 
(second) bistable switch which captures irreversible transitions in hypercortisolism to 
elevated Cortisol levels. 

Conclusions: Model predictive control and inverse bifurcation analysis have been 
previously applied in the simulation-based design of therapeutic strategies for the 
correction of hypocortisolism. Given the HPA model extension presented in this paper, 
these techniques may also be used in the study of hypercortisolism. As an example, we 
show how sparsity enforcing penalization may suggest network interventions that 
allow the return from elevated Cortisol levels back to nominal ones. 
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Background 

The HPA axis is a major part of the neuroendocrine system in mammals and particu- 
larly in humans. A main task of this hormonal network is the regulation of the response 
to physical or mental stress that threatens to disrupt the homeostatic balance of the 
organism. If a stressor is sensed by the nervous system the hypothalamus is stimu- 
lated to produce and secret the corticotropin-releasing hormone (CRH). The secretion 
of CRH causes the anterior pituitary to synthesize adrenocorticotropin (ACTH). ACTH 
then stimulates the adrenal glands to release Cortisol, which down regulates the blood 
concentration of CRH and ACTH via different negative feedback mechanisms [1-4]. 

The HPA axis is subject of intensive research in endocrinology as HPA malfunctions 
are implicated in various pathological conditions. These are often characterized by either 
elevated or insufficient blood Cortisol levels compared to the average healthy human. 
For instance, hypocortisolism (insufficient Cortisol level) is reported in patients suffering 
from the chronic fatigue syndrome and post traumatic stress disorders (cf. [5-9]), whereas 
hypercortisolism (elevated Cortisol level) is observed in depression, dementia or postop- 
erative delirium (cf. [10-15]). Especially in the context of personalized medicine the use of 
modeling and simulation of biological systems for the rational design of treatments and 
drug intervention strategies is more and more recognized [16-20]. For such endeavors 
the integration of biological information of different type into computational, hence pre- 
dictive, models is a prerequisite. The emphasis of earlier HPA modeling approaches with 
ordinary and delay differential equations has been put on self regulatory ultradian and cir- 
cadian oscillatory behavior in [21-27], oscillations in response to an independent external 
pacemaker from the suprachiazmatic nucleus have been described in [28,29]. In compar- 
ison, the article [30] stands out as it incorporates intracellular glucocorticoid receptor 
kinetics which mediate bistable behavior of the HPA axis. Despite its parsimonious char- 
acter the four state ODE model of [30] offers an explanation for hypocortisolism as an 
irreversible biological switch and served in [31] as a basis for the design of a therapeu- 
tic corrections of the HPA dysfunction. In [19] it is shown that the model of [30] also 
exhibits stable limit cycle oscillations, in [32] the four state rate equations of [30] were 
modified in order to fit oscillatory data of patients suffering from post traumatic stress 
disorder. 



A parsimonious HPA model featuring hypocortisolism 

The HPA axis model of [30] captures the basic feedback mechanisms and includes an 
intracellular glucocorticoid receptor GR in the anterior pituitary gland as one of the four 
state variables, see Figure 1. The dynamics of the model are described by the ODE system 
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Figure 1 The HPA axis feedback network. Diagram of the biochemical feedback network of the HPA axis as 
modeled in [30]. Physical or mental stress triggers the secretion of corticotropin releasing hormone CRH in 
the hypothalamus. CRH is transported to the anterior pituitary gland and stimulates the release of the 
adrenocorticotropic hormone ACTH. ACTH signals to the adrenal gland to secrete Cortisol COR. After binding 
to the glucocorticoid receptor GR in the pituitary, Cortisol negatively regulates the production of CRH and 
ACTH. The homo dimerization introduces a positive feedback loop giving rise to bistable behavior in 
accordance with hypocortisolism. 



where the ODE variables [CRH], [ACTH], [GR] and [COR] denote the concentrations 
of CRH, ACTH, GR and Cortisol If GR is bound to Cortisol, the resulting complex is 
denoted by [COR-GR]. As this binding is very fast compared to the other dynamics it is 
assumed to be in equilibrium, i.e. [COR-GR] = [COR] • [GR], The GR dimerization may 
generate an irreversible bistable switch in the model by which the dysfunctional abidance 
of the system at an alternative stable steady state characterized by low Cortisol level can 
be explained. The therapeutic strategy of [31] for the correction of hypocortisolism based 
on (l)-(4) is to force the HPA axis by a further suppression of Cortisol to a point where the 
only stable condition in proximity is the original healthy state. An alternative for support- 
ing the search for pharmaceutical therapies is the manipulation of the feedback network 
by means of inverse bifurcation techniques in order to transform the irreversible switch 
into a reversible one, see [19,33]. Though the model of [30] sheds new light on possi- 
ble causes for hypocortisolism, its dynamic repertoire could not be shown to feature the 
opposite scenario of hypercortisolism. Furthermore, while the model takes the intracel- 
lular glucocorticoid receptor GR into account, it does not consider feedback mechanisms 
known to act on faster time scales. 

Motivation of model extensions 

The above mentioned GR mediated negative feedback on the production of ACTH acts at 
the genomic level. Specifically, Cortisol COR binds to the glucocorticoid receptor GR and 
forms a homodimer (cf. [34,35]). This homodimer is known to work as a transcription 
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factor, enhances the expression of the GR gene and down regulates the expression of the 
ACTH related gene. A consequence of the transcriptional regulation is that the GR feed- 
back control works rather slowly compared to other cellular processes such as a vesicle 
transport, cell signaling or extracellular events such as changes in the CRH blood con- 
centration. For mammalian cells one can expect at least a delay of the down regulation in 
the range of 15 minutes up to 2 hours (c.f. [1]). However, it has been frequently reported 
that Cortisol induced inhibitory effects on anterior pituitary gland cells can already occur 
after a few seconds (c.f. [36,37]), a phenomenon that cannot be explained by means of the 
genomic feedback mechanism. 

The common hypothesis for this observation is that besides the GR receptor a sec- 
ond receptor in the membrane may directly inhibit the release of ACTH when bound to 
Cortisol (cf. [38-42]). In [43] the authors were able to provide evidence of a glucocorti- 
coid receptor (GPCR) in the anterior pituitary cell membrane, possibly enabling the fast 
response related to Cortisol induced inhibition. In the following we incorporate this fast 
and non-genomic feedback mechanism and for this purpose integrate a detailed pituitary 
gland cell model into the HPA network of [30]. As a consequence, the extended dynamic 
repertoire of the resulting HPA axis feedback network also features bistable behavior 
compatible with the dysfunctional state of hypercortisolism and the ability to simulate 
time resolved ACTH stress response data from in vitro AtT-20 pituitary tumor cells. 

Results 

Using mass action and Michaelis Menten enzyme kinetics we developed a mechanis- 
tic ODE system model of the glucocorticoid feedback mechanisms within the anterior 
pituitary gland cell. In particular, we incorporated the glucocorticoid membrane recep- 
tor GPCR found in [43] and distinguished between slow genomic and fast non-genomic 
processes by using different compartments for the nucleus and the cytoplasm. Figure 2 
displays the reaction network graph of the resulting HPA axis model while the details are 
discussed in the Methods section. 

Model comparison with experimental data 

The experimentally observed fast inhibitory effect of Cortisol on the anterior pituitary 
gland cell motivated the model incorporation of the GPCR receptor in order to allow for 
non-genomic, hence fast, signal processing. In order to test the capability of our model to 
capture these fast effects against actual data, we conducted experiments in cultures of the 
murine AtT-20 pituitary cell line (cf. [44]). After administering doses of CRH and Corti- 
sol samples were taken from the medium and analyzed with respect to the total ACTH 
concentration in the cell medium, see Figure 4A for the data obtained for two different 
experimental setups. In the first setup (indicated by the triangle markers in Figure 4A) we 
considered the situation of a single stress stimulus by lOnM CRH at time point zero in 
the absence of Cortisol. We observed an immediate increase of the ACTH secretion, fol- 
lowed by a slight stagnation and then a constantly increasing total ACTH concentration. 
In the second experiment (indicated by the dot markers in Figure 4A) lOnM CRH and 
in addition lOOnM Cortisol were administered. As depicted, this immediately activated 
the fast non-genomic pathway inhibiting the release of ACTH. This led to an increased 
amount of ACTH stored in the vesicles within the cell and thus to a delayed but strong 
release of ACTH after approximately 15 minutes. However, after 30 to 60 minutes the 
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Figure 2 Integration of an anterior pituitary gland cell model into the HPA network. Diagram of the 
biochemical feedback network of the HPA axis after incorporation of the Cortisol mediated feedback 
mechanisms in the anterior pituitary gland cells. The recently found membrane glucocorticoid receptor 
GPCR [43] and the intracellular receptor GR introduced in [30] act on different time scales. The modeling of 
the corresponding feedback mechanisms requires a distinction between slow genomic processes within the 
nucleus and fast non-genomic processes within the cytoplasm. In particular, the dimerized 
Cortisol-glucocorticoid receptor complex (GR complex) and the transcription factors TF associated to the CRH 
receptor CRHR are explicitly considered. See Figure 3 for further details. 






Figure 3 The three core feedback pathways of the (anterior) pituitary gland cell model. Illustration of 
the intracellular pathways of the extended HPA model of Figure 2 showing the three main feedback 
mechanisms (with respect to the involved receptors) in the anterior pituitary gland cells. A) Genomic 
pathway of the GR mediated feedback: Outline of the modeled species involved in the inhibition of POMC 
(Proopiomelanocortin) transcription which corresponds to the inhibition of ACTH production. ACTH is 
obtained via post-translational cleavaging of the POMC protein. The central transcription factor is the 
GR-cortisol complex which down regulates the POMC transcription and enhances the production of the GR 
receptor. B) Non-genomic pathway of the GPCR mediated feedback: The non-genomic negative feedback 
mechanism via the glucocorticoid membrane receptor GPCR, inhibiting the release of ACTH. C) CRHR 
mediated feedback: The positive feedback mechanisms associated with the CRH membrane receptor (CRHR), 
which causes an increased release of ACTH from vesicles and an increased expression of the POMC gene. 
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Figure 4 Experimental and simulated ACTH concentration profiles. The figure shows experimental 
ACTH profiles (black and blue lines) from murine AtT-20 pituitary cell lines described in [44] and simulations 
(green and red lines) after data fitting with different HPA network models. A) The black line corresponds to 
an experiment in which the system was exposed to 10 nM of CRH at time zero, the blue line corresponds to 
an experiment in which - in addition to 1 0 nM - 1 00 nM of Cortisol was administered at time zero. The error 
bars show the 95% confidence intervals for the measurements at the different time points. B) Fit of 
experimental data with the HPA network model of Figure 2. C) Fit of experimental data with the HPA network 
model of Figure 2 after knockout of the membrane glucocorticoid receptor GPCR pathway. D) Fit of 
experimental data with the HPA network model of Figure 1 . The results underpin the need for including the 
membrane glucocorticoid receptor GPCR in models of the pituitary gland cells. 



genomic feedback mechanisms become effective resulting in an overall reduced ACTH 
concentration in the medium. 

For training the model against the experimental data by adaptation of the model param- 
eters we followed a variational regularization approach (cf. [33,45]), see Methods for 
details. Sparsity enforcing penalization was used to eliminate model components less rel- 
evant for reproducing the observed kinetics. Figure 4B shows that the experimental data 
can be fitted by the HPA model of Figure 2 even if the autocatalytic loop for the GR gene 
regulation would be neglected and if the production of the TF (transcription factors for 
CRH) receptors would be considered to be independent of CRH. 

In order to assess the relevance of the featured GPCR receptor for the model we 
knocked out the corresponding pathway by elimination of the respective ODE equations 
from our model and addressed the data fitting task with the reduced model. The result 
in the data space is shown in Figure 4C and indicates that the consideration of the GPCR 
feedback pathway in the model is in fact vital for reproducing the experimentally observed 
dynamics. For comparison we finally run the data fitting procedure with the parsimo- 
nious model (1) - (4) from [30]. The model lacks a GPCR mediated feedback pathway 
for fast response and fails to reproduce the corresponding experimental observation, see 
Figure 4D. 

Exploration of the dynamic repertoire of the model 

One important aspect of [30] is the inclusion of the glucocorticoid receptor GR which 
generates bistability for a model based explanation of hypocortisolism as an irreversible 
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biological switch. As the GR mediated feedback pathway is also a constituent of the 
extended HPA network model of Figure 2, the capability of our model to reflect the 
bistable behavior of (l)-(4) is an imperative requirement. Tools [46,47] for studying the 
dynamic repertoire of ODE systems with respect to bistability or oscillations solve opti- 
mization problems that involve the eigenvalues of the Jacobian matrix f x (x,q) (with/ 
denoting the right hand side of an ODE model) and search for regions in the parame- 
ter space that give rise to limit-point or Hopf bifurcations. Figure 5B displays a resulting 
bifurcation diagram for the extended HPA model from Figure 2 with the stress level of 
Equation (8) as the bifurcation parameter. It gives an explanation of the dysfunctional 
alteration from a steady state of normal Cortisol level (upper branch of the equilibrium 
curve) to a second steady state with low Cortisol level (lower branch of the equilibrium 
curve) in response to a temporary stressor. If the abscissa of the limit point LP2 lies below 
the basal stress level, say for instance 0.1, the return to the upper branch is hampered and 
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Figure 5 Dynamic Repertoire of the HPA network model. The figure displays different types of qualitative 
dynamic behavior featured by the HPA axis model from Figure 2 which is larger than for other HPA network 
models reported earlier in the literature. The corresponding model parameter settings have been found 
using inverse problems and bifurcation discovery tools. The bifurcation diagrams (B) and (D) show the 
steady state concentration of Cortisol in dependence on the stress level. For a certain stress region the system 
exhibits two stable steady states (indicated by solid lines) and one unstable steady state (dashed line). 
A) Simulation of a dysfunctional stress response of Cortisol in hypocortisolism. Though the stress is only 
temporary (red period), the system is driven to and eventually caught in an alternative stable steady state 
characterized by an insufficient Cortisol level. B) Irreversible bistable switch explaining the behavior shown in 
(A). Once the system is driven from the upper steady state branch it is trapped at the lower branch due to a 
limit point LP 2 abscissa smaller than the basal stress level 0.1 . C) Simulation of a dysfunctional stress response 
of Cortisol in hypercortisolism. Though the stress is only temporary (red period), the system is driven to and 
eventually caught in an alternative stable steady state characterized by an excessive Cortisol level. 
D) Irreversible bistable switch explaining the behavior shown in (C). Note the inverted alignment of Cortisol 
steady state branches and alternative location of limit points in comparison to (B). E) Simulation of limit cycle 
oscillations in the hormone concentrations of CRH, ACTH and Cortisol (relative to the maximum Cortisol level). 
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the HPA axis is trapped at low Cortisol as observed in hypocortisolism. Figure 4A shows 
a corresponding Cortisol concentration profile obtained by simulation of the model from 
Figure 2. 

Further explorations of the dynamic repertoire of the extended HPA model revealed a 
region in the parameter space which allows for an alternative bistable behavior that can 
be linked to the pathological condition of hypercortisolism. Figure 5D provides an exam- 
ple of a corresponding bifurcation diagram in which the qualitative scenario of Figure 5B 
is inverted. If a stressor drives the HPA axis from the lower Cortisol branch to the upper 
one, it abides at elevated Cortisol levels even after removal of the stressor if the abscissa of 
the limit point LP2 lies below the basal stress level, again say for instance 0.1. Figure 5C 
shows a corresponding Cortisol concentration profile obtained by simulation of the model 
depicted in Figure 2. While the positive feedback mediated by the glucocorticoid receptor 
GR dimerization has been already suggested in [30] as the cause for the bistable behavior 
associated to hypocortisolism, our in-silico pathway knock out tests on the extended HPA 
model indicate that the positive feedback on the genomic level mediated by the CRH hor- 
mone is vital for the bistability associated to hypercortisolism. Without this component, 
we could not find a parameter region giving rise to a bifurcation diagram as displayed in 
Figure 5D. Similarly, the results of our numerical analysis suggest that the dynamic reper- 
toire of the parsimonious model of [30] excludes the bifurcation behavior of Figure 5D, 
possibly due to its simplistic description of the CRH pathway. 

Finally, we explored the parameter space of the extended HPA network model with 
objective functions for the detection of Hopf bifurcations and found parameter regions 
allowing for stable limit cycle oscillations, see Figure 5E. The displayed temporal sequence 
in which the core species CRH, ACTH and Cortisol reach their maximum value is in 
accordance with the established cascade in release of the HPA hormones. While oscilla- 
tory behavior for the parsimonious model was obtained in [30] in response to repeated 
stress pulses, limit cycle oscillation for (l)-(4) were found in [19] and for a modified ver- 
sion also in [32] . However, to our knowledge the HPA model shown in Figure 2 is the first 
one able to capture two distinct types of bistability (as reflected in Figures 5B and 5D) as 
well as autonomous oscillatory behavior. 

Model based correction of Hypercortisolism 

Based on the parsimonious model of Figure 1 from [30] a therapeutic strategy for the 
correction of hypocortisolism has been suggested in [31]. This therapy aims at perturbing 
the HPA axis trapped at the low Cortisol branch by a further suppression of Cortisol in 
order to enable the return to a healthy Cortisol level. An alternative idea for the treatment 
of hypocortisolism has been presented in [33] which aims at an intervention with the 
reaction network in order to turn the dysfunctional irreversible switch into a reversible 
one by shifting the critical limit point in the bifurcation diagram. As the extended HPA 
model of Figure 2 also captures a switching scenario associated to hypercortisolism, the 
inverse bifurcation approach of [33,48,48,49] can also be applied to study which parts of 
the HPA network need to be targeted in order to allow the return from the abnormal 
hypercortisol steady state back to a healthy Cortisol level. 

For the purpose of illustration we suppose that a representative q° out of the parameter 
region compatible with the switching behavior of Figure 5D has been chosen (for instance 
by fitting of corresponding patient data) in which the abscissa S2(q°) of the limit point 
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LP2(q°) lies below the basal stress level S = 0.1 (on a fictive scale). Then, a short but 
strong stress signal may drive the system from a normal Cortisol level to an abnormal 
elevated one where it abides even if the stress drops back to basal level, see Figure 5C. A 
therapeutic qualitative change from an irreversible to a reversible steady state transition 
can be achieved by interfering with the HPA model such that the abscissa of the limit 
point LP2(q°) is shifted to a value with > 0.1 while, e.g., maintaining the abscissa 
value of the second limit point LP\ (q°), see Figure 5 A. 

A corresponding interference strategy can be derived by solving the nonlinear inverse 
problem 



where F maps the correction x onto the abscissas S\(q° + x) and S2(q Q + x) of the limit 
points LP\ (q° + x) and LP2 (q° + x) defined by the HPA model with altered parameter 
vector q° + x. In order to keep the number of necessary interventions with the HPA net- 
work as low as possible a sparse solution of (5) is sought for by solving the optimization 
problem 



with (p as in (35) and the sparsity promoting choice p < 1. 

Figure 6C displays the solution of (6) (in relation to q°) giving rise to the re-engineered 
bifurcation diagram of Figure 6A (green curve) and allowing the return to normal Corti- 
sol levels, see Figure 6B. The parameters V3 and d-j to be altered concern the feedback 
controls related to the hormone CRH. A decrease of V3 and an increase of d-j correspond 
to a reduction of the sensitivity of the cells with respect to extracellular CRH. Moreover, 
an increase of the dissociation/Michaelis-Menten constant K$ reduces the sensitivity of 
the TFs transcription factor regulation by the CRH-CRHR related pathways. Further- 
more, the increase of /c t i2 calls for an increase of the translation rate of GR mRNA in 
order to further diminish the effect of extracellular CRH due to higher availability of GR. 
The intervention strategy shown in Figure 6C also suggests to alter the parameters 
ds and ^2 in order to reduce the production of CRH and to increase the sensitivity with 
respect to the feedback via Cortisol. As these parameters concern processes outside of the 
anterior pituitary gland cells, an extension of HPA model by a mechanistic description of 
the hypothalamus is necessary in order to enable a more detailed interpretation of those 
parameter changes. 

Conclusions 

Motivated by the experimental detection of a glucocorticoid membrane receptor GPCR 
we have built a mathematical model of the anterior pituitary cell that differentiates 
between genomic and GPCR mediated non-genomic signaling pathways. After the inte- 
gration of our model into a literature framework of the hypothalamus pituitary adrenal 
axis and the training against stress induced ACTH concentration profiles of murine 
AtT-20 pituitary tumor cells we obtained a dynamic repertoire capturing two types of 
bistable switches as well as stable limit cycle oscillations. While the first bistable switch 
has been previously associated to the diseased state of hypocortisolism and utilized in 
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Figure 6 Targeting the irreversible switch of the HPA axis model for hypercortisolism. A) The blue 
equilibrium curve corresponds to the irreversible switch of hypercortisolism depicted also in Figure 4B. A shift 
of the limit point LP 2 (q°) to the region right to the basal stress level turns the irreversible switch into a 
reversible one (green equilibrium curve). This design goal can be formulated and solved using inverse 
problems techniques. B) Comparison of the dysfunctional stress response with the normal stress response 
after network manipulation allowing the return to normal Cortisol levels following a short stress period (red 
area). C) Display of the necessary model parameter corrections x (relative to the nominal values ofq°) in 
order to achieve the transition from an irreversible to a reversible switch as shown in Figure 5A. Due to the 
use of a sparsity promiting penalty function only a few parameters need to be altered. 



the development of therapeutic strategies, to our knowledge the second bistable switch is 
novel in the context of HPA axis modelling. Along similar lines it can be linked to another 
dysfunctionality of the HPA axis referred to as hypercortisolism in which a continued 
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abidance of the axis at elevated Cortisol levels is observed even after complete removal of 
the stressing factor. Our in-silico knockout studies in particular indicate that the positive 
feedback on the genomic level mediated by the CRH hormone is responsible for the gen- 
eration of the hypercortisolism switch. Using inverse bifurcation analysis with sparsity 
enforcing penalization we finally detected key model parameters that need to be altered 
in order to turn the irreversible switch into a reversible one then featuring a healthy stress 
response characterized by a return to normal Cortisol levels. 

While systems biology strongly advertises the use of modeling and simulation in the 
design of future drugs and therapies, one current challenge is to appropriately picture the 
complexity of the biological system at hand in the defining mathematical equations and 
their parameters. Though with respect to HPA axis modeling we reached a certain level of 
detail by introducing both genomic and non-genomic signaling pathways in the pituitary 
cell, we only adopted the parsimonious description of the link between the pituitary and 
the hypothalamus as well as the adrenal from the literature. Furthermore, the information 
content of the available ACTH concentration did not allow to sufficiently constrain the 
model parameters such that we so far have to content ourselves with parameter space 
regions corresponding to different qualitative model behaviour. Nevertheless, even at its 
current stage our model offers to endocronologists a different view on hypercortisolism 
as an irreversible bistable switch while the outlined mathematical methodologies will also 
apply to more realistic models and improved data sets. Finally, the model may also serve 
as a tool for the design of a model predictive control based treatment of hypercortisolism 
in which the HPA axis is driven by a - though at first glance counter-intuitive - additional 
administration of Cortisol to a point where the only stable condition in proximity is the 
original healthy state, compare with [31]. 

Methods 

Model development 

In order to distinguish between slow genomic and fast non-genomic processes we use the 
modeling concept of compartments and separate the nucleus from the cytoplasm of the 
pituitary gland cell. Moreover, we base our model on the principles of Michaelis-Menten 
enzyme kinetics ([50-52]) in order to deal with physical meaningful parameters. 

We refrain from a detailed description of gene regulation, RNA processing or cleaving 
and rather concentrate on the rate limiting processes such as transport/diffusion in and 
out of the nucleus. We apply this modeling strategy not only to the glucocorticoid medi- 
ated feedbacks but also to the stimulating pathways involving the hypothalamic hormone 
CRH. CRH does not enter the pituitary gland cells but interacts with a specific mem- 
brane receptor (CRHR) only. Still, also the binding of CRH induces both a genomic and 
a non-genomic feedback (cf. [53-55]). As this positive feedback via CRH counteracts the 
negative feedback via Cortisol on the genomic as well as on the non-genomic level, also 
the catalyzing mechanisms related to CRH have to be modeled in more detail than in [30] 
for the purpose of differentiating between genomic and non-genomic effects in the ante- 
rior pituitary gland cells. Concerning the molecular mechanisms of the anterior pituitary 
cells, our model involves three different pathways which are associated with the three 
main receptors, namely the intracellular glucocorticoid receptor GR and the two mem- 
brane receptors (GPCR and CRHR) with respect to Cortisol and CRH, see Figure 2 and 
Figure 3. 
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Extracellular species and the controls related to the adrenals and the hypothalamus 

Regarding the feedback controls related to the adrenals and the hypothalamus we follow 
the approach of [30] in order to insure consistency and comparability as far as possible. 
In particular, we consider the following ODE describing the extracellular concentration 
of Cortisol, with an additional term reflecting the diffusion (diffusion constant /cdiffi) of 
Cortisol into the anterior pituitary cells: 

d[COR ex ] V in vi [ACTH ex ] 

-77 = 77 ^diffl ([COR in ] -[COR ex ] ) + ~ ; r 7 ____ ; - d\ [COR ex ] . (7) 

at V ex Ai + IAClhLexJ 

Here V ex and V- m are the volumes of the extracellular compartment and the intracellu- 
lar compartment. Analogously to [30] we model the temporal development of the CRH 
concentration by 

d[CRH] = k fa + stress) 

dt l + [COR ex ]//C3 

The GR mediated feedback control 

One significant extension over [30] in our modeling approach concerns the regulatory 
pathway related to the intracellular glucocorticoid receptor GR and in particular to the 
transcription factor, being the homodimer of the GR-COR complex, see Figure 3A. The 
following equation describes the intracellular concentration of Cortisol 

d[C ° Rin] = fcdiffl ([CORe X ] -[COR in ] ) - <k [COR in ] 

dt (9) 

+ F ([(GR-COR) 2 ,nu] , [CORin] , [GR] ) . 

The change in the amount of Cortisol in the cytoplasm (CORin) is subject to the 
transport of the dimerized GR-COR complex into the nucleus. This is described by 

Wurrcrn^ 1 rrnrc 1 \rm\ ■ 7, / *2« ^grc ^Rdim [ (GR-COR) 2 ,n U ] 
F ([ ( GR-COR) 2 ,n U ] , [CORin] - [GR] ) .= k ^ [CORta p. [GR] +[co ^ HGR]2 



fr 2in [CORin] -[GR] 
([CORin] +[GR] ) 



(10) 



which follows from the assumption of the fast binding of Cortisol and GR and the quick 
formation of the homodimer. A detailed derivation is provided in the Additional file 1. 

The dimerized GR-COR complex is transported to the nucleus, where it acts as a 
transcription factor and binds to the respective DNA sites. In particular it binds in the 
promoter region of the GR gene and enhances the transcription of the GR RNA sequences 

d[pmGR] v 8 [(GR-COR) 2 , n „] 

— — = [pmGR] +V7 + /C 9+ [(GR-COR) 2 , nu ] - dl5 [pmGR] ' (U) 

which then are transported to the cytoplasm and translated to the GR protein. 

The second gene regulated by the GR-COR complex is the Proopiomelanocortin gene 
POMC, which yields the ACTH peptide after post-translational cleavage. While the 
POMC gene is up-regulated by certain transcription factors (TFs) related to the CRH 
stimulus, it is down-regulated by the GR-COR complex. It has been observed that these 
two transcription factors counteract via competitive inhibition (cf. [53,56]), i.e. the bind- 
ing of the GR-COR complex prevents the catalyzation of the POMC expression and thus 
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down regulates the transcription. This can be described by the following equation for the 
transcribed RNA 

d[pmPOMC] v 6 [TFs nu ] 

= v 5 + 



dt K 7 (l+[ (GR-COR) 2 , nu ] /K 8 ) +[TFs nu ] 

- /ctrsi [pmPOMC] -d u [pmPOMC] . (12) 
The GPCR mediated feedback control 

The second pathway we consider in our model is the negative feedback on the release of 
ACTH via the glucocorticoid membrane receptor (GPCR), see Figure 3B. Concerning the 
production of the GPCR receptor we assume no regulation by means of other compounds 
and hence disregard the processes on the genomic level With respect to the complex 
formation (COR-GPCR) we suppose a (quasi) equilibrium as it has been experimentally 
observed that two ligands bind to the receptor and exhibit positive cooperativity (cf. [43]). 

The GPCR mediated feedback control only regulates the release of ACTH from intra- 
cellular vesicles. Again this mechanism counteracts the positive stimulus by CRH which 
is also detected via a membrane bound receptor (CRHR). It is well known that the ante- 
rior pituitary gland cells show a fast response to exposure with Cortisol and CRH and 
that this is achieved by signaling cascades, directly relating the stimuli to the fusing of 
the vesicles (cf. [36,37]). Consequently, we assume that the release of ACTH directly 
depends on the amount of formed receptor ligand complexes, where the catalysis of the 
CRH-receptor complex and the inhibitory effect of the cortisol-receptor complex are 
modeled as competitive inhibition. The following ODE describes the extracellular ACTH 
concentration 

d[ACTH ex ] = Vin friex [CRHR-CRH].[ACTH in ] 

dt V ex 7<4 (l+[GPCR-(COR) 2 ] /K 5 ) + CRHR-CRH] 

(13) 

The CRHR mediated feedback control 

The last pathway we consider concerns the feedback mechanism related to CRH and the 
respective membrane receptor CRHR, see Figure 3C. The main effects are the stimu- 
lated release of ACTH and the enhanced expression of the POMC gene. The genomic 
effects of the CRH mediated feedback involve several transcription factors (the Nur 
factors, Tpit/Tbxl9, Pitxl, cf. [53]), which catalyze the expression of the POMC gene. 
However, for keeping the complexity moderate we only incorporate a single CRH related 
transcription pathway into our model and denote the representative factor by TFs, see 
Equation (12). The alternative would be to model each of the transcription pathways sep- 
arately at the cost of additional equations and parameters to be identified, for details 
on transcription modeling at different levels see, e.g., [57]. Consequently, we also model 
the catalyzed production of TFs only in a coarse manner, where we account for the fast 
signaling of the binding of CRH onto the production of TFs by a hill factor, leading to 



d[TFs in ] V nu v 4 [CRHR-CRH]^ 

= — /Oex [TFs nu ] -fan [TFs in ] +^ + [ ^^ " *2 [TFs in ] . 

(14) 
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Concerning the fast response to CRH stimulation it is the common hypothesis that 
the release of ACTH is regulated via a signaling cascade which is induced by the bind- 
ing of CRH to the CRHR receptor in the cell membrane. Signaling cascades are very fast 
compared to other cellular processes. In fact the response to Cortisol or CRH has been 
observed already after a few micro seconds (cf. [40,54,58,59]). Thus we assume an imme- 
diate effect of the membrane receptor complex formation on the release of ACTH, as 
modeled in the ODE (13). 

Figure 2 gives a graphical outline of the complete model which is mathematically 
described by means of the following non-linear ODE system with 46 modeling parameters 
in 15 equations 



d[COR ex ] Vin, nrr%1> 1 rr _ D n vi [ACTH ex ] 

j: = — fcdiffi ([COR in ] - [CORex])+ ^ , -di [COR ex ] 

at V ex AiH-[AClH ex J 

(15) 
(16) 



d[CRH] k s (k s b + stress) 



dt 

d [ACTHex] 

dt 

d[COR in ] 
dt 

d[ACTH in ] 
dt 



d[GPCR] 
dt 

d[CRHR] 
dt 

d[GR] 
dt 



d[(GR-COR) 2 , nu ] 
dt 



dt 

d[mGR] 
dt 



1+ [COR ex ] /Kz 



- d 2 [CRH] 



k lex [CRHR-CRH] -[ACTH in ] 



Vex /<4 (l+[GPCR-(COR) 2 ] /K s ) +[CRHR-CRH] 
- d 3 [ACTH e x] (17) 



= ^diffl ([CORex] — [CORin] ) - d 4 [CORin] 

+ F([ (GR-COR) 2 , nu ] , [CORm] , [GR] ) 

= km [mPOMC] -d 5 [ACTHi n ] 

Ariex [CRHR-CRH] -[ACTH in ] 
~ /<4 (l+[GPCR-(COR) 2 ] /K 5 ) +[CRHR-CRH] 

= V2-ek [GPCR] 
= v 3 -d 7 [CRHR] 



(18) 

(19) 

(20) 
(21) 



= k t n [mGR] +F ([(GR-COR) 2 , nu ] , [COR in ] , [GR] ) - d s [GR] 



(22) 



k 2 in [(GR-COR) 2)in ]-A: 2 ex [(GR-COR) 2 , nu ]-d 9 [(GR-COR) 2 , nu ] 

(23) 



d[mPOMC] Vn, 



Vin 
Vnu 

V in 



*trsi [pmPOMC] -d w [mPOMC] 
^trs2 [pmGR] — dn [mGR] 



(24) 
(25) 
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d [ 1 1 s in ] - nu ^ [ T p Snu ] _/ C3 . n [TFs in ] -d u [TFs in ] 



dt V [r 



d[TFs n 



dt 

d[pmPOMC] 
dt 



v 4 [CRHR-CRH]^ , 

+ —it —r 26 

/C^+tCRHR-CRH]^ 

: /r 3in [TFs in ] -feex [TFs nu ] -rfi3 [TFs nu ] (27) 



Jttrsi [pmPOMC] +v 5 - rfi4 [pmPOMC] 

, V 6 [TFs nu ] (28) 

K 7 (l+[(GR-COR) 2 , nu ]//<C 8 )+[TFs nu ] 

d[pmGR] | v 8 [(GR-COR) 2> nu] , r _ D1 _ 

= " [pmGR] + v 7 + 7 , 9+[(GR , CQR)2 nu] " *s [pmGR] , (29) 

where 

W UrK m m i rrnR i r r pi\. 7, A 2ex ^GRC *GRdim [ (GR-COR) 2 , nu ] 
F ([ (GR-COR) 2 , nu] , [COR in ] , [GR] ) .= * [CORin p . [GR] +[CORjn] . [GR]2 

*2ln[COR ta ]-[GR]\ 



([CORi„] +[GR] ) 



with A: = 



Model parameters 

We now comment on the role of the parameters in our ODE model We have three sub- 
sidiary parameters describing the average volumes of the three compartments V m , V ex 
and V nu . These parameters are phenomenological and are fixed in the subsequent 
analysis. 

Another class of parameters are the constants related to diffusion and transport of the 
compounds. In our model /<diffi> kiex> fein> I<2ex> I<3ex and /<3i n comprise the list of these 
parameters. We do not explicitly model active transport, as this would involve additional 
regulatory species and mechanisms and consequently would lead to an even more com- 
plex model. Subsequently, we assume that the rate of diffusion/transport is constant on 
average. 

The largest portion of model parameters is related to degradation processes whose 
rates are assumed to be proportional to the concentration of the respective species. The 
associated rate constants are denoted by d\, d<i, . . . , <ii4 and d\^. 

As a consequence of the considered enzymatic reactions and ligand-receptor complexes 
the model includes several dissociation constants. Moreover, we have some constants 
arising due to the Michaelis-Menten approximations, i.e. the maximum reaction rates 
and the Michaelis-Menten constants. They are closely related to the dissociation con- 
stant of the respective complex and approximate them if the rate of the product formation 
dominates the degradation of the enzyme-substrate complex. /<GRdim> ^GC2> #crc and 
feRC denote the dissociation constants of the considered ligand-receptor complexes. K\, 
7<9 denote the Michaelis-Menten constants of the respective enzyme-substrate 
complexes, vi, k s , k s /c s b and V4 denote the corresponding limiting rates of the Michaelis- 
Menten approximations. 
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Our description of the genomic mechanisms involves rates for transcription (V5, V6, V7 
and v$), translation (k t \i and ktn) and transport (/ctrsi and /<trs2) of mRNA. We neglect 
the post-processing of mRNA as our focus is on the respective timescales and as the 
time delay is mainly due to the translocation. The constants V2 and V3 denote the pro- 
duction rate of the two membrane receptors. As we have not considered their regulatory 
mechanism we lumped the transcription and translation rate together. 

Model training against data 

The central idea of variational regularization is to augment the least-squares term that 
measures the mismatch between the data and the model with a penalizing term (p that 
guarantees the continuous dependence of the parameter solution on the data 



\\x(q,pi) - y[\\ 2 + \\x(q,p 2 ) - y S 2 L +a<p(q) -> min (31) 

(q,Pl,P2) 

subject to: x(q, Pk) = f(x, q) (32) 

x(q t p k ) = p k ke{l,2} (33) 

c(p k )<0 ke{l,2] (34) 



The ODE system (32) represents the HPA axis model of Figure 2 along with its param- 
eter vector q of length 46. p\ and p 2 in (33) denote the initial conditions corresponding 
to the two different experimental settings as described above and the two resulting data 
sets y\ and y\ displayed in Figure 5A. While the initial concentration of some species are 
known from the experimental setup and can be taken into account by the constraints (34), 
the remaining components of pi, p 2 have to be estimated from the data as well. As the 
available amount of data is not sufficient to uniquely determine the model parameters, 
we have chosen a two step approach for finding the model components most relevant for 
reproducing the data. First, we applied classical Tikhonov regularization [45] and solved 
(31) with p = 2 in the penalty term 

i 

using a combination of global and local optimization routines. The resulting parameter 
vector q tlk then served as an initial guess for solving (31) with a j?-value less than one. 
Sparsity enforcing i p functionals with the choice 0 < p < 1 aim at solutions whose 
less important components are driven to zero [60-62] and have been used in the context 
of parameter identification in [63,64] . For nonlinear inverse problems, the regularization 
properties of (35) with 0 < p < 1 have been analyzed in [65,66]. Figure 7 displays the 
parameter vector cf obtained for p = 0.9 in relation to q tlk . The near zero components 
of (f correspond to reactions and pathways that can rather be neglected in fitting the 
experimental data, see Figure 5B. In particular, the autocatalytic loop for the GR gene 
regulation, represented by the parameters k tvs , V7, vg, V9, d\5, d\\, seems negligible with 
respect to the fast dynamics. Furthermore, the near zero value of K& indicates that the 
production of the TFs receptors is independent of CRH, i.e. the genomic part of the CRH 
related feedback. 
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v ° v model parameters 

Figure 7 Comparison of parameter solutions for data fitting. Plot of q*/q tlk (in componentwise sense) 
with the parameter vectors q* and q tlk as inferred from the experimental data of Figure 4A in a two step 
approach. The result shows that certain pathways of the model are of minor relevance for the reproduction 
of the experimentally observed dynamics as given in Figure 4B. 



Origin of experimental data 

Murine anterior pituitary AtT-20 cells (passage 12-24) were seeded and maintained as 
previously described [44]. After 92 hours of cell growth AtT-20 cells were exposed to 
doses of lOnM CRH and in addition lOOnM Cortisol. The supernatant was carefully 
removed from the cell layer after 1 minute to 1 hour of administration and cen- 
trifuged for 10 minutes with 800xg at 37° C. The concentration of ACTH-molecules 
in the cell-free supernatant was detected using Fluorescence Correlation Spectroscopy 
Immunoassay [44]. 

The AtT-20 cells (ATCC no. CCL-89) were purchased from the American Type Culture 
Collection (ATCC, Manassas, USA). Cell culture media and reagents were from Sigma 
(Sigma- Aldrich Inc., St. Louis, MI). 
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Additional file 1: Detailed derivation of the mathematical model. 
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